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Abstract 

Generalized  polynomial  chaos  (gPC)  methods  have  been  successfully  applied 
to  various  stochastic  problems  in  many  physical  and  engineering  fields.  However, 
realistic  representation  of  stochastic  inputs  associated  with  various  sources  of  un¬ 
certainty  often  leads  to  high  dimensional  representations  that  are  computationally 
prohibitive  for  classic  gPC  methods.  Additionally  in  the  classic  gPC  methods, 
the  gPC  bases  are  determined  based  on  the  probabilistic  distribution  of  stochastic 
inputs.  However,  the  stochastic  outputs  may  not  share  the  same  probabilistic  distri¬ 
bution  as  the  stochastic  inputs.  Hence,  the  gPC  bases  may  not  be  the  optimal  bases 
for  such  systems,  which  causes  the  slow  convergence  of  gPC  methods  for  such 
stochastic  problems.  Here  we  present  a  general  framework  that  integrates  the  adap¬ 
tive  ANOVA  decomposition  technique  and  the  data-driven  stochastic  method  to 
alleviate  both  of  the  two  limitations.  To  handle  high-dimensional  stochastic  prob¬ 
lems,  we  investigate  the  use  of  adaptive  ANOVA  decomposition  in  the  stochastic 
space  as  an  effective  dimension-reduction  technique  for  high-dimensional  stochas¬ 
tic  problems.  Three  different  ANOVA  adaptive  criteria  are  discussed. 

To  improve  the  slow  convergence  of  gPC  methods,  we  use  the  data-driven 
stochastic  method  (DDSM)  which  was  developed  by  Cheng-Hou-Yan  in  [5].  This 
method  has  an  offline  computation  and  an  online  computation.  In  the  offline  com¬ 
putation,  optimal  gPC  bases  are  obtained  by  Karhunen-Loeve  (K-L)  expansion 
of  the  covariance  matrix  of  stochastic  outputs  obtained  by  ANOVA-based  sparse- 
grid  PCM.  In  the  online  computation,  a  Galerkin-projection  based  gPC  method 
with  the  optimal  bases  developed  in  the  offline  computation  is  employed,  which 
greatly  speeds  up  the  convergence.  Numerical  examples  are  presented  for  one- 
,  two-dimensional  elliptic  PDE  with  random  coefficients,  and  a  two-dimensional 
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1  Introduction 

There  has  been  growing  interest  and  significant  progress  over  the  past  decades  on  mod¬ 
eling  complex  physical,  and  engineering  systems  with  uncertainties.  However,  in  the 
simulation  of  real  physical  randomly  heterogeneous  systems,  a  realistic  representation 
of  stochastic  inputs  associated  with  various  sources  of  uncertainty  often  leads  to  high 
dimensional  representations  that  are  computationally  prohibitive  for  many  numerical 
methods,  such  as  generalized  polynomial  chaos  (gPC)  methods  [38]  and  probabilis¬ 
tic  collocation  method  (PCM)  [37].  In  this  paper,  we  present  a  general  framework 
that  combines  the  ANOVA  decomposition  technique  with  the  recently  developed  data- 
driven  stochastic  method  (DDSM)  [5]  to  alleviate  the  curse  of  dimensionality  and  slow 
convergence  issues  of  gPC  methods.  For  stochastic  problems  with  high  stochastic  di¬ 
mensions,  we  employ  the  functional  ANOVA  (ANalysis-Of-VAriance)  method  [15,  2] 
as  a  dimension-reduction  technique.  The  ANOVA  decomposition  was  introduced  by 
Fisher  [8].  Later,  Hoeffding  in  1948  successfully  apply  ANOVA  decomposition  to 
study  U-statistics  [16].  This  method  is  motivated  by  the  observation  that  for  many  real 
physical  systems,  only  a  relatively  small  number  of  stochastic  dimensions  are  impor¬ 
tant  and  will  significantly  impact  the  outputs  of  the  stochastic  systems.  ANOVA  has 
also  been  used  for  uncertainty  quantification  in  [36]  and  was  employed  in  gPC  for  solv¬ 
ing  high-dimensional  stochastic  PDF  systems  in  [9,  24,  25,  4,  11,  39].  In  [9]  ANOVA 
was  integrated  with  a  multi-element  PCM.  In  [24],  an  adaptive  version  of  ANOVA  is 
developed  to  automatically  detect  the  important  dimensions.  In  [39],  adaptive  ANOVA 
methods  based  on  three  different  adaptive  criteria  were  proposed  and  compared. 

The  ANOVA  decomposition  results  in  a  set  of  low-dimensional  sub-problems  in 
stochastic  space.  A  sparse-grid  PCM  is  used  to  solve  these  low-dimensional  sub¬ 
problems  efficiently.  The  PCM  was  first  introduced  by  Tatang  and  McRae  [33] .  Rrecently 
Xiu  and  Hesthaven  [37]  have  used  a  Lagrange  polynomial  interpolation  to  construct 
high-order  stochastic  collocation  methods.  The  properties  of  PCM  were  extensively 
studied  in  the  past  10  years.  In  [26,  27,  1],  the  errors  of  integrating  or  interpolating 
functions  with  Sobolev  regularity  were  analyzed  for  Smolyak  constructions  based  on 
the  one-dimensional  nested  Clenshaw-Curtis  rules.  In  [27],  the  degree  of  exactness 
of  the  Smolyak  quadrature  using  the  Clenshaw-Curtis  and  Gaussian  one-dimensional 
rules  was  investigated.  In  [37],  the  efficiency  of  the  Clenshaw-Curtis-based  sparse  grid 
stochastic  collocation  was  demonstrated  by  comparing  it  with  other  stochastic  methods 
on  an  elliptic  problem.  In  2003,  Gerstner  and  Griebel  [12]  introduced  the  dimension- 
adaptive  tensor  product  quadrature  method.  In  [10],  sparse  grid  collocation  schemes 
were  applied  to  solve  stochastic  natural  convection  problems.  In  [21,  22,  18,  13],  a 
multi-element  PCM  was  employed  to  study  the  random  roughness  problem,  stochas¬ 
tic  compressible  flow  and  plasma  flow  problems.  In  [14,  23],  an  adaptive  hierarchical 
sparse  grid  collocation  algorithm  has  been  developed. 
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In  the  classic  gPC  methods,  the  gPC  bases  are  determined  based  on  the  probabilistic 
distribution  of  stochastic  inputs.  DDSM  is  inspired  by  the  fact  that  the  stochastic  out¬ 
puts  may  not  share  the  same  probabilistic  distribution  as  the  stochastic  inputs.  Hence, 
the  distribution  used  for  sparse-grid  PCM  may  not  be  the  optimal  distribution  to  repre¬ 
sent  the  solutions  of  such  systems,  which  causes  the  slow  convergence  of  sparse-grid 
PCM  for  such  stochastic  problems.  To  overcome  this  difficulty,  we  use  DDSM  [5]  to 
obtain  a  set  of  problem-dependent  optimal  gPC  bases  to  greatly  speed  up  the  conver¬ 
gence.  The  number  of  gPC  modes  to  achieve  a  specified  accuracy  used  in  the  ANOVA- 
based  DDSM  is  much  smaller  than  the  classic  gPC  method,  which  greatly  reduces  the 
computational  cost.  In  this  work,  an  ANOVA-based  DDSM  is  developed,  which  can  be 
considered  as  a  stochastic  extension  of  the  Proper  Orthogonal  Decomposition  (POD) 
methods  [31,  35],  the  original  DDSM  of  Cheng-Hou-Yan  [5],  and  the  Multiscale  Finite 
Element  methods  [17].  DDSM  has  an  offline  computation  and  an  online  computation. 
In  the  offline  computation,  optimal  gPC  bases  are  obtained  by  Karhunen-Loeve  (K-L) 
expansion  of  the  covariance  matrix  of  stochastic  outputs  obtained  by  ANOVA-based 
sparse-grid  PCM.  In  the  online  computation,  a  Galerkin-projection  based  gPC  method 
with  optimal  bases  developed  in  the  offline  computation  is  employed,  which  accelerates 
the  convergence  rate  considerably. 

The  remainder  of  the  paper  is  organized  in  the  following  way;  section  2  describes 
the  standard  ANOVA  decomposition,  section  3  presents  the  DDSM,  Section  4  intro¬ 
duces  the  ANOVA-DDSM.  Section  5  shows  the  results  obtained  using  ANOVA-DDSM 
and  section  6  summarizes  the  findings. 


2  Sparse-Grid  Based  Probabilistic  Collocation  Method 


The  general  procedure  for  the  PCM  approach  is  similar  to  MC  simulations,  with  a 
difference  in  selecting  the  sampling  points  and  corresponding  weights.  The  procedure 
consists  of  three  main  steps: 

1.  Generate  collocation  points  in  probability  space  of  random  parameters  as 
independent  random  inputs  based  on  a  quadrature  formula; 

2.  Solve  a  deterministic  problem  at  each  collocation  point; 

3.  Estimate  the  solution  statistics  using  the  corresponding  quadrature  rule. 


{u{x,t))  =  (  u{x,t,^)p{^)d^ 


(1) 


Jr 

Nc 


(2) 


a{u){x,  t) 


(3) 


(4) 
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where  p(^)  is  the  probabilistic  distribution  function  (PDF)  of  random  variable  Nc  is 
the  number  of  quadrature  points,  {^fc}  is  the  set  of  quadrature  points  and  {wfc}  is  the 
corresponding  set  of  weights,  which  are  the  combination  of  quadrature  weights  in  each 
random  dimension.  In  the  second  step  of  the  PCM  approach,  as  for  MC,  any  existing 
code  can  be  used  to  solve  deterministic  flow-and-transport  equations.  Extensive  re¬ 
views  on  the  construction  of  quadrature  formulas  may  be  found  in  [6]  and  [7].  Below, 
we  provide  a  brief  review  of  two  different  methods  for  selecting  collocation  points. 

In  this  work,  we  use  the  Smolyak  formula  [32]  to  construct  the  collocation  point 
set,  which  is  a  linear  combination  of  tensor  product  formulas,  and  the  resulting  point  set 
has  a  significantly  smaller  number  of  points  than  the  full  tensor  product  set.  Recently, 
researchers  [37,  19,  20]  have  used  Lagrange  polynomial  interpolation  to  construct  high 
order  stochastic  collocation  methods  based  on  sparse  grids  using  the  Smolyak  formula 
[32].  Such  sparse  grids  do  not  depend  as  strongly  on  the  dimensionality  of  the  random 
space  and  as  such  are  more  suitable  for  applications  with  high-dimensional  random 
inputs.  Detailed  description  on  how  to  build  the  collocation  point  set  can  be  found  in 
[37,  19,  20]. 

3  Data-Driven  Stochastic  Method 

Consider  the  stochastic  PDE 


£{x,uj)u{x,uj)  =  f{x) 


(5) 


where  x  G  D,uj  G  and  C{x,u!)  is  a  stochastic  differential  operator.  The  stochastic 
ingredient  resides  in  the  differential  operator  C{x,  to)  while  f{x)  is  purely  determinis¬ 
tic.  The  authors  heuristically  argue  that  the  bases  only  depend  on  £(x,  w).  Generally 
speaking,  the  bases  can  be  obtained  by  solving  (5)  for  one  particular  f{x)  in  offline 
part.  Then  this  set  of  bases  could  be  used  to  solve  other  problems  with  different  com¬ 
plicated  f{x). 

Inspired  by  multiscale  finite  element  method  [17]  and  proper  orthogonal  decom¬ 
position  (POD)  method  [31,  35],  Cheng-Hou-Yan  proposed  an  new  algorithm,  called 
Data-Driven  Stochastic  Method  [5].  The  authors  tried  to  build  up  gPC  bases  under 
which  the  stochastic  solutions  have  a  sparse  decomposition  based  on  Karnunen-Loeve 
(K-L)  expansion.  Eor  m(x,  w)  G  L'^{D  x  O),  its  K-L  expansion  reads  as  follows 


(6) 


where  {A^}  and  {(/)i(x)}  are  the  eigenpairs  of  the  covariance  kernel  C'(x,  y),  i.e. 


(7) 


and  {^i(a;)}  are  random  variables  defined  as 


V  Ai  J  D 


/  (u(x,  w)  —  ]E[u])0i(x)(ix, 

Jd 


(8) 


4 


Moreover,  these  mutually  uncorrelated  random  variables  {■Ci(w)}  satisfy 


(9) 


In  numerical  practice,  only  a  finite  series  expansion  is  adopted,  depending  on  the  decay 
rate  of  the  eigenvalues. 


N 


(10) 


For  a  given  covariance  function,  the  decay  rate  of  the  eigenvalues  depends  inversely 
on  the  correlation  length.  Long  correlation  length  implies  that  the  random  process  is 
strongly  correlated  and  results  in  a  fast  decay  of  the  eigenvalues.  A  weakly  correlated 
process  has  short  correlation  length  and  results  in  a  slow  decay  of  the  eigenvalues. 
Under  some  assumptions  [30],  the  eigenvalues  in  K-L  expansion  decay  exponentially 
(or  sub-exponentially)  fast  in  dimension  d  =  1  (or  d  >  1). 

This  algorithm  consists  of  offline  part  and  online  part.  In  offline  part,  a  set  of 
gPC  bases  {Ai{uj)}  are  obtained  based  on  K-L  expansion.  In  online  part,  the  set  of 
gPC  bases  are  used  to  solve  a  set  of  coupled  deterministic  PDFs  of  coefficients  in  the 
stochastic  expansion.  It  should  be  noticed  that  in  this  framework  the  set  of  gPC  bases 
are  problem-dependent. 

3.1  Offline  computation 

The  purpose  of  offline  computation  is  to  obtain  a  complete  set  of  gPC  bases  based  on  K- 
L  expansion  of  the  stochastic  solutions.  Instead  of  using  Monte  Carlo  (MC)  simulation, 
which  only  has  a  convergence  rate  of  1/ \/iV,  we  combine  sparse  grid  method  and  K-L 
expansion  in  order  to  expedite  the  computation  in  high  dimensional  problem.  In  this 
part,  we  obtain  not  only  the  gPC  bases  numerically,  but  also  a  variety  of  statistical 
information  on  this  set  of  bases  {Ai{u!)},  such  as  expectation  of  a{'x,uj)Ai{u!)Aj{uj) 
and  moments 


In  order  to  construct  the  covariance  function  C{x,  y),  we  need  to  sample  stochastic 
solution  u(x,  w).  First,  we  expand  m(x,  uj)  using  polynomial  chaos  expansion: 


(11) 


Then  by  orthogonality  of  {'0j(u;)},  we  have 


Uj{x)  =  /  u{x,uj)tjjj{uj)di 

Jn 


(12) 


The  integral  (12)  can  be  approximated  by  sparse  grid  method. 


(13) 
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where  {uJi}  are  Smolyak  sparse  grid  points,  {£ti(x,  Wi)}  are  corresponding  solutions, 
and  {wi}  are  integration  weights  associated  with  {uii}.  Once  {■Uj(x)}  are  obtained, 
we  can  sample  N  solutions  m(x,  w)  using  K-L  expansion  (11).  Then  the  mean  and 
covariance  are  computed 


^  i^l 
1  ^ 

c'lx,  y)  =  X]  « W“(y)- 

i=l 

Afterward,  the  first  M  eigen-pairs  are  solved 


(14) 

(15) 


X^(j)i{x)=  /  C{x,y)cj),{y)dy,  i  =  (16) 

JD 

Finally,  the  gPC  bases  {Ai{uj)}  are  obtained  by 

Ai{u})  =  ^=l  {u{x,uj)  -  u{x))(j)i{x)dx.  (17) 

V  Xi  J  D 

Remark  1  It  is  easy  to  verify  that  each  basis  Ai{uj)  has  mean  zero  and  {Ai(u;)}^Q 
are  mutually  orthogonal. 


3.2  Online  computation 

In  this  part,  we  only  solve  deterministic  equations  since  all  the  statistical  information 
has  been  obtained  in  the  offline  part.  This  means  that  the  online  computation  could  be 
fast.  The  bases  Ai{u!)  spans  a  finite-dimensional  subspace  in  L^(0).  Therefore,  we 
can  project  the  stochastic  solution  it(x,  ui)  onto  this  subspace,  i.e. 

M 

m(x,w)  Ri  (18) 

i=0 

For  simplicity  of  notation,  Aq  =  1  and  Uo(x)  =  E[u(x,  w)]. 

In  order  to  obtain  the  coupled  deterministic  equations,  Galerkin  projection  is  uti¬ 
lized.  Multiplying  (5)  by  Aj(uj)  and  taking  expectation  on  both  sides  give  us 

M 

'^E[C{yi,uj)Ai{uj)Aj{uj)]ui{x)  =E[f{x)Aj{u})],  j  =  (19) 

i^O 

4  ANOVA  Expansion 

In  statistics,  ANOVA  method  can  be  used  to  describe  the  interactions  between  a  large 
number  of  variables  while  only  few  samples  are  available.  The  same  idea  can  be 
adopted  in  the  interpolation  and  integration  of  high  dimensional  problems  as  well  as 
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stochastic  systems.  For  most  well-defined  physical  system,  only  relatively  low  order 
correlations  of  the  input  variables  are  expected  to  be  important  for  the  output  of  the  sys¬ 
tem.  The  ANOVA  expansion  utilizes  this  property  and  at  each  new  level  of  ANOVA 
expansion,  higher  order  correlation  effects  of  the  input  variables  are  accounted  for. 
Consider  a  Lebesgue  integrable  multivariate  stochastic  function  /(Y)  :  — >■  M  and  d 
is  the  dimension  of  stochastic  space  we  are  interested  in.  ANOVA  expansion  represents 
/(Y)  as  finite  hierarchical  correlated  functions  of  input  variables  in  the  form  of 

d 

/(Y)  =  /o  +  ^  ^  (20) 

s=i  ii<"-<is 

or  equivalently 

/(Y)  =  /o+  ^  fMJ+  E  fn.J2iYn^Y,J  +  --- 

+  /i,2,..,d(Yi,Y2,---  ,Yd).  (21) 

We  call  fj^  {Yj^ )  the  first  order  term,  fj^.  {Yj^, ,  Y,, )  the  second  order  term,  etc.  . 

The  ANOVA  components  have  the  following  properties: 

1.  The  constant  term  is  the  mean  of  function,  that  is 

fo=  [  f{Y)df,{Y),  (22) 

which  means  that  all  higher  order  components  have  mean  zero 

/  /j,,...  ,jJn{Y)  =  0  for  1  <  s  <  d.  (23) 

Jr** 

2.  The  other  important  property  of  ANOVA  expansion  is  the  orthogonality  among 
its  terms 

/  fn,-,jJk„-MMY)=0,  (24) 

if  (jii  ■  ■  ■  I  Js)  ^  (fci,  •  •  •  ,  ki).  This  is  the  direct  consequence  of  (23). 

3.  The  variance  of  /  is  the  sum  of  variance  of  all  component  functions 

=E  E (25) 

s=l  |s|=s 

It  is  worth  pointing  out  that  equation  (25)  holds  only  when  the  measure  used  in  the 
calculation  of  variance,  i.e.  the  integral  with  Lebesgue  measure,  is  the  same  as  that  in 
ANOVA  decomposition. 


7 


Remark  2  It  could  be  extremely  expensive  to  compute  AN OVA  decomposition  for  high 
dimensional  /(Y).  Therefore  Dirac  measure  is  adopted  instead  of  Lebesgue  measure, 
i.e.  dp,{'Y)  =  i5(Y  —  c)dY,  c  £  The  special  point  c  is  termed  anchor  point.  How¬ 
ever,  it  is  difficult  to  calculate  anchor  point  c  such  that  /g  =  /(c)  =  /(Y).  In  this 
paper,  we  take  anchor  point  c  to  be  the  mean  of  random  variable  Y  as  an  approxima¬ 
tion.  In  this  case,  the  property  (22)  and  (23)  do  not  hold  any  more.  Additional  terms  of 
ANOVA  decomposition  are  needed  to  improve  the  mean. 

The  measure  dp,{Y)  determines  the  particular  form  of  each  component  function 
following  the  notation  in  [29,  28].  We  introduce  a  projection  operator  T^g  :  F'^  — ^  fI®I 

iP./(Yg)  :=  [  f{Y)dpi\,{Y)  (26) 

where  :=  dpi{Yi). 

Therefore,  each  term  fs  can  be  recursively  defined  by 

/s(Yg)  =  iPs/(Yg)  -  ^  /t(Yt).  (27) 

tCs 


4.1  Adaptive  ANOVA 


When  the  nominal  dimension  of  the  stochastic  problem  increases,  the  computational 
complexity  of  the  standard  ANOVA  becomes  prohibitive  to  evaluate  all  the  terms.  For 
example,  for  nominal  dimension  N  =  100,  the  number  of  terms  for  2nd  order  ANOVA 


decomposition  needed  to  calculate  is  1  + 


=  5051.  Nevertheless, 


in  many  stochastic  problems,  most  of  the  interactions  among  different  dimensions  are 
usually  weak  and  have  little  contribution  to  the  stochastic  outputs.  This  means  that 
the  active  dimension  of  those  stochastic  problems  is  small.  Therefore,  an  adaptive 
approach  can  be  employ  to  solve  those  problems  efficiently  without  losing  much  accu¬ 
racy. 

There  are  many  “adaptive”  approaches  and  the  one  we  employ  in  this  paper  is 
obtained  by  replacing  the  nominal  dimension  by  an  active  dimension,  i.e.,  we  modify 
(21)  to  be 


/(Y)«/o+  ^  MYn)+  E  + 

ji<Di  (lij'2)  6.7^2 

+  E  (28) 

In  practice,  Di  is  usually  set  to  be  N  and  v  =  2.  That  is  all  first  order  terms  are 
calculated  so  that  some  criterion  could  determine  the  set  J'2  according  to  the  property 
of  the  specific  problem.  However,  the  criteria  in  [25,  39]  does  not  specify  how  many 
the  active  dimensions  there  are  a  priori.  Hence  there  may  be  redundant  dimensions 
which  are  actually  not  active  dimensions.  Motivated  by  the  DDSM  method,  the  fast 
decay  of  eigenvalues  in  (7)  can  serve  as  an  appropriate  indicator  of  the  convergence  of 
K-L  expansion  in  terms  of  optimal  gPC  bases. 


First,  we  describe  two  popular  adaptive  criterion  listed  in  [25,  39].  Then  a  new 
criterion  is  proposed  in  criterion  3. 

Criteria  1:  Let  Ti  =  which  is  the  sum  of  the  variance  of  all  the 

first-order  terms.  Assume  that  the  first  order  terms  are  sorted  such  that  (J^{fj)  is  mono- 
tonically  decreasing.  The  active  dimension  D2  should  satisfy: 

£>2 

(29) 

i=i 


where  p  is  a  proportionality  constant  with  0  <  p  <  1  and  very  close  to  1.  This  criterion 
is  similar  to  the  criterion  used  in  [3]  where  cr^(/)  instead  of  Ti  is  used  on  the  right 
hand  side  of  (29)  and  p  is  set  to  be  0.99.  The  set  can  be  found  by  computing 


'^31,32 


31,32) 


(30) 


and  bounding  with  a  predefined  error  threshold  6*2. 

Criteria  2;  Ma  and  Zabaras  use  the  mean  of  component  function  fj  as  the  indicator 
to  decide  the  active  ANOVA  terms  [25].  Let 


E(/,) 

/o  ’ 


(31) 


where  the  predefined  error  threshold  9i  is  used  to  bound  rjj,  i.e.,  rjj  <  61  for  some  J. 
If  r]j  are  monotonically  decreasing  with  respect  to  j  and  we  set  D2  >  J,  then  (31)  can 
equivalently  be  written  as 


£>2 


N 


^E(/,)>p^E(/,), 

i=i 


(32) 


where  p  =  1  —  9i.  Then,  for  a  further  selection  of  the  second-order  terms,  Ma  and 
Zabaras  also  use  the  mean  of  the  component  functions: 


9 31, 3 2 


^^(/jl,j2) 

EfiomY 


(33) 


where  rij^  j^  is  bounded  by  a  predefined  error  threshold  02- 

Criteria  3:  Combining  all  above,  we  propose  a  new  criterion,  which  truly  selects 
active  dimensions.  First,  eigenvalues  Ai  of  covariance  kernel  C{3i,y)  in  equation  (7) 
are  obtained  and  sorted,  i.e.  Ai  >  A2  >  •  •  • .  Then  we  pick  the  number  of  active 
dimensions  D2  such  that 


Ai/Ai  >  A,  i  —  1,  •  •  •  ,1)2,  (34) 

where  A  is  a  predefined  threshold.  The  dimensions  with  largest  D2  either  variance  or 
mean  are  selected  as  active  dimensions.  Further  selection  of  the  second  dimensions 
J^2  could  be  accomplished  by  (30)  or  (33)  accordingly.  Our  numerical  examples  of  Id 
and  2d  Elliptic  PDE  with  random  coefficients  show  that  this  criterion  is  much  more 
efficient  by  the  fact  that  D2  is  much  less  than  that  in  Criteria  1  or  2. 
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Remark  3  When  we  employ  the  above  criteria  to  applications  we  replace  the  mean 
and  deviation  of  component  function  fj  with  their  L2  norm  values  on  the  physical 
domain. 


5  (Adaptive)  ANOVA-based  Data-Driven  Stochastic  Method 
(ANOVA-DDSM) 


In  this  paper,  we  solve  high  dimensional  stochastic  problems  by  taking  advantage  of 
(adaptive)  ANOVA  method  and  data-driven  stochastic  method.  For  large  dimension 
d  ^  1,  the  rate  of  convergence  in  many  stochastic  methods,  such  as  Monte  Carlo 
method,  Wiener  Chaos  Expansion  method  (WCE),  etc.  ,  deteriorates  drastically.  This 
poses  a  numerical  challenge  because  it  requires  a  huge  number  of  simulations  of  the 
underlying  deterministic  system.  This  is  the  well-known  curse  of  dimensionality.  The 
data  driven  stochastic  method  also  has  this  problem  since  the  number  of  gPC  bases 
M  increases  fast  if  the  dimension  of  the  stochastic  problem  is  large.  The  following 
algorithm  is  proposed  to  deal  with  high  dimensional  stochastic  problems  efficiently. 

Algorithm  1  ((Adaptive)  ANOVA-DDSM  algorithm) 

1.  Expand  stochastic  solution  u(x,  Y)  in  ANOVA  decomposition 


i(x,Y)=  Mo(x) -f  ^  Ujj(x,YjJ -t-  +  ••■ 

ii=i  i<ii<i2<ci 


(35) 


or  adaptive  ANOVA  decomposition 

Di 


t(x,Y)Ri  uo(x)  +  ^  -f  +  , 

ii=i  (iij2)e.7^2 

(36) 


where  the  set  T2  is  selected  according  to  Criterion  3. 

2.  Solve  for  the  mean  term  uo(x)  which  satisfies  a  deterministic  equation  by  re¬ 
placing  random  variables  Y  with  anchor  point  c  in  the  stochastic  PDE, 

/:(x,c)uo(x)  = /(a;).  (37) 


3.  Eor  each  high  order  term  j^{Yj^,Yj^,  -  ■  ■  DDSM  is  utilized  to 

calculate  the  solution  efficiently  with  different  deterministic  forcing  term  f{x). 
Eor  different  f{x),  the  same  gPC  bases,  which  have  been  constructed  in  the 
offline  part  for  each  term  in  ANOVA  decomposition,  can  be  used  repeatedly. 

Remark  4  In  practice,  only  second  order  or  third  order  terms  are  needed  in  ( adaptive) 
ANOVA  expansion  for  good  accuracy  of  mean  and  variance. 

Remark  5  In  calculating  each  term  in  ANOVA  expansion,  several  gPC  bases  could  be 
enough  for  accuracy  requirement.  This  essentially  expedites  our  computation. 
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6  Numerical  examples 


As  explained  in  previous  sections,  we  can  reduce  a  high  dimensional  stochastic  prob¬ 
lem  into  a  system  of  low  dimensional  problems  using  (adaptive)  ANOVA  decomposi¬ 
tion.  For  each  low  dimensional  problem,  it  is  significantly  efficient  to  solve  each  low 
dimensional  problem  in  terms  of  CPU  time  and  memory  cost. 

6.1  ID  Stochastic  Elliptic  PDE  in  4d  Random  Space 

We  test  the  ANOVA-DDSM  algorithm  on  a  ID  stochastic  Elliptic  PDE  with  4  sources 
of  randomness  as  follows 


(38) 


u{x,  w)  =  0,  X  =  0  and  x  =  1. 


The  stochastic  coefficient  a(x,  w  in  Eq.  (38)  is  chosen  to  be 


4 


(39) 


where 


C  =  [0.1,  0.12,  0.2,  0.15],  D  =  [1.2,  2.3,  3.1,  4.3] 


(40) 


and  {^i}  are  uniform  i.i.d.  random  variables  in  [0,1].  We  run  10^  realizations  of 
Monte  Carlo  as  the  exact  solution  umc(x,  w).  The  DDSM  gPC  bases  {Ai{uj)}f£i  are 
constructed  using  /(x,  y)  =  I  and  then  are  used  to  solve  the  stochastic  Elliptic  PDE 
with  /(x,  y)  =  sin(27rx)  +  5  sin(47r2/).  The  relative  L2  error  of  mean  is  defined  as 


\\u{x)  -  Umc{^)\\l2 
\\umc{^)\\l2 


(41) 


and  relative  L2  error  of  variance  is 


||Var(u)  -  Var(uMc)|lL2 

||Var(uMc)||  L2 


(42) 


The  relative  error  of  mean  is  shown  in  Eig.  1 .  The  constant  term  does  not  repre¬ 
sent  the  mean  of  the  stochastic  solution  well  since  the  anchor  point  c  is  not  optimal. 
Therefore,  additional  terms  are  necessitated  to  improve  the  accuracy  of  ANOVA  ap¬ 
proach.  Note  that  the  stochastic  coefficient  a(x,  uj)  is  an  additive  function  of  random 
variables  Then  the  first  order  expansion  is  enough  to  represent  the  solution  and 

its  mean  [25],  which  is  shown  in  Eig.  1.  This  plot  indicates  that  increasing  the  ex¬ 
pansion  order  does  not  improve  accuracy  dramatically.  Moreover,  it  only  has  limited 
effect  on  the  accuracy  of  mean  by  including  more  gPC  bases  (or  modes)  in  expansion 
(18).  Eirst  order  ANOVA  expansion  and  4  or  5  stochastic  modes  are  good  enough  to 
approximate  the  mean  of  solution.  However,  higher  order  ANOVA  decomposition  and 
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Figure  1 ;  Relative  error  of  mean  of  one  dimensional  elliptic  PDF  with  4  random  vari¬ 
ables.  The  exact  solution  is  computed  with  10^  realizations. 


more  stochastic  modes  are  required  to  represent  the  variance  well.  As  shown  in  Fig. 
2,  second  order  of  ANOVA  decomposition  and  6  bases  are  needed  to  approximate  the 
variance  of  stochastic  solution.  Contrast  to  the  mean,  increasing  the  order  of  ANOVA 
decomposition  and/or  number  of  stochastic  bases  do  improve  the  accuracy  of  variance. 
In  practice,  unless  otherwise  stated,  we  take  second  order  ANOVA  expansion  and  6 
gPC  bases  in  our  calculations. 


6.2  ID  Stochastic  Elliptic  PDE  in  High  Dimensional  Prohahilistic 
Space 

In  this  subsection,  we  consider  ID  stochastic  Elliptic  PDE  (38)  in  100  dimensional 
probabilistic  space.  The  stochastic  coefficient  a(x,uj)  now  reads  as 

too 

a(x,  uj)  =  '^  +  1),  (43) 

i=l 

where  are  uniform  random  variables  in  [0, 1]  and  Ci  G  (0, 0.001)  and  Di  G  (0, 10) 
are  randomly  generated.  We  run  10®  realizations  of  Monte  Carlo  as  the  reference  solu¬ 
tion  umc{^,  w).  The  DDSM  gPC  bases  {Aiiu))}^^  are  constructed  using  f{x,  y)  =  1 
and  then  are  used  to  solve  the  stochastic  Elliptic  PDE  with  f{x,y)  =  sin(27ra;)  + 
5  sin(47ry). 
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Figure  2:  Relative  error  of  variance  of  one  dimensional  elliptic  PDF  with  4  random 
variables.  The  exact  solution  is  computed  with  10^  realizations. 


Table  1 :  Comparison  of  adaptive  ANOVA-DDSM  using  different  Criteria 


relative  error  of  mean 

relative  error  of  cr 

#  of  active  terms 

Criterion  1 

(p  =  0.9,  ?7j,j2=1E-5) 

2.04E-2 

2.14E-2 

151 

Criterion  3 

A  =  lE-3 

A  =  lE-5 

A  =lE-3 

A  =lE-5 

A  =lE-3 

A  =lE-5 

3.60E-2 

2.23E-2 

2.21E-2 

2.15E-2 

111 

146 

The  decay  of  A^/Ai  is  shown  in  Figure  3.  Only  5  dimensions  are  active  if  A  =1E- 
3  and  10  dimensions  if  A  =lE-5.  However,  if  we  pick  p  =  0.9,  the  numbers  of 
active  dimensions  are  54  and  87  for  Criterion  1  and  2  respectively.  Thus,  the  numbers 
of  second  order  terms  needed  to  be  computed  are  1431  and  3741  respectively.  The 
computations  of  all  these  terms  are  necessary,  although  the  number  of  active  terms 
in  second  order  could  be  further  shrinked  by  imposing  some  threshold  rjj^  for  the 
relative  variance  or  mean  in  Criterion  1  or  2.  In  contrast,  in  Criterion  3,  the  active 
dimensions  have  been  identified  beforehand  for  some  predefined  threshold  A.  Thus  the 
number  of  second  order  terms  needed  to  compute  is  significant  less. 

The  comparison  using  Criteria  1  and  3  is  listed  in  Table  1.  From  this  table,  it  is 
clear  that,  using  different  criteria,  the  numbers  of  final  active  terms  are  almost  the  same 
and  so  are  the  relative  errors  of  mean  and  standard  deviation  a. 

The  profile  of  mean  computed  by  second  order  adaptive  ANOVA-DDSM  using 
Criterion  3  is  compared  with  MC  method  in  Fig.  4.  It  can  be  seen  that  the  profile 
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Figure  3:  The  decay  of  sorted  relative  eigenvalues  of  covariance  kernel  for  a  ID  Elliptic 
PDE  with  random  coefficient  in  100  probabilistic  dimension. 


matches  the  reference  solution  very  well.  However,  the  standard  deviation  of  ANOVA- 
DDSM  is  a  little  lower  than  MC  as  shown  in  Eig.  5.  Higher  order  adaptive  ANOVA 
decomposition  can  definitely  improve  the  accuracy  of  standard  deviation.  But  it  brings 
considerably  more  amount  of  computations. 


6.3  Horn  Problem 


In  this  subsection,  we  extend  our  method  to  2-d  Helmholtz  equation  in  random  media, 
the  planar  acoustics  horn  problem  described  in  detail  in  [34].  The  full  domain  is  de¬ 
picted  in  Eig.  6,  which  comprises  both  the  horn  proper  and  a  large  circular  segment. 
The  governing  equations  for  the  (complex)  pressure  are  then 


VMx,y,w)  +  k'^{l+n'^{x,y,u}))p{x,y,uj)  =  0, 

with  boundary  conditions 


dp 

dn 


ikp 

dp 

dn 


p{x,y,u}) 


0  on  r  1 , 

0  onr2, 

f{x,y)  on  Eg, 


(44) 
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Figure  4:  Comparison  of  mean  of  one  dimensional  elliptic  PDF  with  100  random 
variables  by  Adaptive  ANOVA-DDSM  using  Criterion  3  and  Monte  Carlo  simulation. 
A  =lE-5  and  10  active  dimensions. 


Figure  5:  Comparison  of  standard  deviation  of  one  dimensional  elliptic  PDF  with  100 
random  variables  by  Adaptive  ANOVA-DDSM  using  Criterion  3  and  Monte  Carlo  sim¬ 
ulation.  A  =lF-5  and  10  active  dimensions. 
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where  n  is  the  unit  outer-pointing  normal  of  the  boundary,  k  is  the  wave  number  and 
n^(x,  y,  oj)  is  the  random  reflectivity  of  the  media.  In  this  example,  the  random  reflec¬ 
tivity  of  the  media  is  chosen  to  be 


4 


where  {'Ci(w)}  are  i.i.d.  uniformly  distributed  random  variable  in  [0, 1]  and  functions 
{tpiix,  y)}  are  given  by 


tpi{x,  y)  =  sit?{2'kx  +  9i)  sin^(27ry  +  62) 
ip2{x,  y)  =  sin^(47ra;  -f  63)  sin^(47ry  -f  64) 
l/’a  I  y)  =  sin^  (67ra;  +  9^ )  sin^  (47ry  +  9q  ) 
ip4(x,  y)  =  sin^(67ra;  +  9^)  sin^(67ry  +  9s) 


with  phase  {9i}  being  randomly  generated. 

The  DDSM  gPC  bases  {Ai{uj)}f^i  are  constructed  using  f{x,y)  =  1  and  then 
are  used  to  solve  the  stochastic  horn  problem  with  f{x,  y)  =  sin(27ra;)  sin(27rj/)  +  2. 
In  this  example,  k  is  taken  to  be  0.7.  Only  6  gPC  bases  and  second  order  ANOVA 
decomposition  are  used.  In  Fig.  7,  the  comparison  of  real  part  of  mean  contours  of 
pressure  by  ANOVA-DDSM  and  MC  methods  is  shown.  As  expected,  the  mean  could 
be  approximated  accurately  even  using  only  second  order  ANOVA  decomposition  in 
this  example.  However,  the  difference  of  standard  deviation  is  noticeable,  especially 
outside  the  horn  proper  and  the  region  right  of  center.  Again,  this  is  mainly  due  to  the 
lower  order  of  ANOVA  decomposition.  Accuracy  of  deviation  could  be  compensated 
by  employing  higher  order  ANOVA  decomposition  and/or  more  gPC  bases.  Adaptive 
ANOVA  method  could  be  utilized  if  computational  cost  is  highly  concerned  and  high 
accuracy  is  required. 

6.4  2D  Stochastic  Elliptic  PDE  in  High  Dimensional  Probahilistic 
Space 

Finally,  We  consider  a  2D  stochastic  Elliptic  PDE  with  a  random  coefficient  in  50 
dimensional  probabilistic  space 


-V  •  {a{x,y,uj)\7u{x,y,uj))  =  f{x,y), 
u{x,y,^)\aD  =  0. 

The  stochastic  coefficient  a{x,  y,  to)  is  defined  as 

50 


(45) 


(46) 


where  are  uniform  random  variables  in  [0, 1]  and  the  constant  parameters  are  ran¬ 
domly  generated  Di  G  [0,  0.01],  Ei,  Fi  G  [5, 10].  Adaptive  ANOVA-DDSM  methods 
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Figure  6;  Geometry  of  horn  problem 


Figure  7:  Comparison  of  real  part  of  mean  contour  in  horn  problem.  Left:  ANOVA- 
DDSM  method.  Right:  Monte-Carlo  simulation. 
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Figure  8:  Error  of  real  part  of  mean  contour  in  horn  problem. 


Figure  9:  Comparison  of  standard  deviation  contour  in  horn  problem.  Left:  ANOVA- 
DDSM  method  .  Right:  Monte-Carlo  simulation. 
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Figure  10;  Error  of  standard  deviation  contour  in  horn  problem. 


19 


Table  2:  Comparison  of  adaptive  ANOVA-DDSM  using  different  Criteria 


relative  error  of  mean 

relative  error  of  a 

#  of  active  terms 

Criterion  1 

(p  =  0.99,  ?7jij2=1E-5) 

3.60E-3 

5.70E-2 

201 

Criterion  3 

A  =  lE-3 

A  =  lE-4 

A  =lE-3 

A  =lE-4 

A  =lE-3 

A  =lE-4 

4.60E-3 

4.40E-3 

6.72E-2 

6.18E-2 

61 

96 

using  different  criteria  are  employed  to  solve  this  problem.  The  DDSM  gPC  bases 
{Ai{uj)}f^i  are  constructed  using  f{x,  y)  =  I  and  then  are  used  to  solve  the  stochas¬ 
tic  Elliptic  PDE  with  /(x,  y)  =  sin(27ra;)  sin(27ry)  +  2.  Monte-Carlo  results  with  10® 
realizations  are  served  as  reference  solutions.  The  decay  plot  of  eigen-values  of  covari¬ 
ance  kernel  in  Eig.  1 1  indicates  that  there  are  5  active  dimensions  if  Xi/Xi  >lE-3  and 
10  active  terms  if  Xi/Xi  >lE-4  in  Criterion  3.  Eor  Criterion  1,  the  parameter  p  is  set 
to  be  0.99,  which  results  in  25  active  dimensions.  Then  total  300  second  order  terms 
are  computed.  It  requires  201  terms  in  total  for  =1E-5  to  be  calculated  towards 

approximating  stochastic  outputs.  In  contrast,  to  achieve  the  same  order  of  accuracy 
for  both  mean  and  standard  deviation,  there  are  only  61  and  96  terms  needed  in  total 
for  A  =lE-3  and  A  =lE-4,  respectively.  The  comparison  results  are  summarized  in 
Table  2. 

In  Eig.  12,  the  contour  of  mean  computed  by  adaptive  ANOVA-DDSM  using  Cri¬ 
terion  3  with  10  active  dimensions  is  plotted  and  compared  with  that  by  MC.  In  terms 
of  mean,  adaptive  ANOVA-DDSM  performs  well.  To  further  demonstrate  the  conver¬ 
gence  of  mean,  we  plot  the  error  contour  of  mean  in  Eig.  13.  This  shows  clearly  that 
the  magnitude  of  error  is  of  order  10“^.  In  addition,  we  plot  the  mean  contour  lines 
of  stochastic  solutions  given  by  adaptive  ANOVA-DDSM  and  MC  in  Eig.  14.  Again, 
the  contours  match  very  well.  All  the  above  have  shown  that  ANOVA-DDSM  captures 
mean  accurately. 

In  Eig.  15,  the  contour  of  standard  deviation  computed  by  adaptive  ANOVA- 
DDSM  using  Criterion  3  with  10  active  dimensions  is  plotted  and  compared  with  that 
by  MC.  Adaptive  ANOVA-DDSM  can  capture  the  correct  pattern  of  the  standard  devi¬ 
ation.  However,  the  magnitude  is  slightly  lower.  This  is  indicated  in  Eig.  16,  where  the 
error  of  standard  deviation  is  plotted.  In  order  to  further  illustrate  this,  the  contour  lines 
of  standard  deviation  are  shown  in  Eig.  17.  The  closer  to  the  center  of  the  domain,  the 
larger  the  solution  is.  In  this  plot,  the  contour  lines  of  adaptive  ANOVA-DDSM  are 
closer  to  the  center  than  those  of  MC,  which  means  that  the  deviation  calculated  by 
adaptive  ANOVA-DDSM  is  smaller.  This  is  mainly  due  to  the  following  two  facts. 
Eirst,  only  10  active  dimensions  in  adaptive  ANOVA  decomposition  are  considered 
and  this  is  just  an  approximation  of  the  standard  ANOVA  decomposition.  Secondly 
and  more  importantly,  only  second  order  decomposition  is  employed.  Using  higher  or¬ 
der  adaptive  ANOVA  decomposition  could  improve  the  accuracy  of  standard  deviation, 
but  it  would  take  much  more  computational  cost. 
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Figure  11:  The  decay  of  sorted  relative  eigenvalues  of  covariance  kernel  for  a  2D 
Elliptic  PDE  with  random  coefficient  in  50  probabilistic  dimension. 
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Figure  12:  Comparison  of  mean  contour  of  2D  stochastic  Elliptic  PDE.  Left:  adaptive 
ANOVA-DDSM  method  using  Criterion  3.  Right:Monte-Carlo  simulation.  A  =lE-4 
and  10  active  dimensions. 


mean  error 


Eigure  13:  Mean  error  contour  of  2D  stochastic  Elliptic  PDE. 
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Figure  14:  Contour  comparison  of  mean  contour  of  2D  stochastic  Elliptic  PDE. 


Eigure  15:  Comparison  of  standard  deviation  contour  of  2D  stochastic  Elliptic  PDE. 
Left:  adaptive  ANOVA-DDSM  method  using  Criterion  3.  Right:  Monte-Carlo  method. 
A  =lE-4  and  10  active  dimensions. 
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Figure  16:  Standard  deviation  error  contour  of  2D  stochastic  Elliptic  PDE. 


Eigure  17:  Contour  comparison  of  standard  deviation  contour  of  2D  stochastic  Elliptic 
PDE. 
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7  Conclusion  and  Discussion 


In  this  paper,  a  novel  adaptive  ANOVA-based  data-driven  method  is  developed  for 
solving  high-dimensional  stochastic  elliptic  equation  arising  from  various  applications, 
such  as  the  randomly  heterogeneous  porous  media  flow  problem.  The  developed  method 
has  an  offline  computation  and  an  online  computation.  In  the  offline  computation, 
adaptive  ANOVA  decomposition  technique  is  applied  to  adaptively  decompose  the 
original  high  dimensional  problem  into  a  set  of  low-dimensional  sub-problems.  By 
modeling  the  behavior  of  stochastic  systems  with  only  the  first  few  lower-order  terms 
of  the  high-dimensional  input,  adaptive  ANOVA  is  able  to  efficiently  represent  the 
output  response  to  the  high-dimensional  inputs  with  specified  good  accuracy.  This 
results  in  a  set  of  low-dimensional  sub-problems  in  stochastic  space,  which  are  effi¬ 
ciently  solved  by  sparse-grid  PCM.  Numerical  examples  have  shown  that  solving  the 
set  of  low-dimensional  sub-problems  is  more  efficient  than  solving  the  original  prob¬ 
lem.  Three  different  ANOVA  adaptive  criterion  are  discussed.  Numerical  tests  indicate 
that  the  third  adaptive  criteria  gives  the  best  approximation  with  minimal  computational 
cost. 

Numerical  examples  involving  both  one-,  two-dimensional  elliptic  PDE  with  ran¬ 
dom  coefficients,  and  a  two-dimensional  Helmhotz  equation  in  random  media  (Horn 
problem)  have  been  conducted  to  verify  the  accuracy  and  efficiency  of  the  developed 
adaptive  ANOVA-based  DDSM  method.  In  the  offline  computation,  for  stochastic 
problems  with  fixed  number  of  stochastic  dimension,  the  number  of  component  func¬ 
tions  needed  in  adaptive  ANOVA  decomposition  depends  on  the  important  dimensions 
with  respect  to  the  stochastic  outputs  and  the  variance  of  the  stochastic  inputs.  For 
real  physical  high-dimensional  stochastic  problems  with  up  to  500  —  600  stochastic 
dimensions  [9,  25],  adaptive  ANOVA  decomposition  integrated  with  sparse-grid  PCM 
can  achieve  much  better  convergence  rate  than  both  the  Monte  Carlo  and  sparse-grid 
PCM.  However,  it  is  worthwhile  to  note  that  adaptive  ANOVA  decomposition  may  not 
be  recommended  to  approximate  mathematical  functions  where  all  dimensions  are  im¬ 
portant.  Additionally,  sparse-grid  PCM  is  determined  based  on  the  probabilistic  distri¬ 
bution  of  stochastic  inputs.  However,  due  to  the  nonlinearity  of  the  complex  stochastic 
systems,  the  numerical  solutions  may  not  share  the  same  probabilistic  distribution  as 
the  stochastic  inputs.  Hence,  the  distribution  used  for  sparse-grid  PCM  may  not  be  the 
optimal  distribution  to  represent  the  solutions  of  such  systems,  which  causes  the  slow 
convergence  of  sparse-grid  PCM  for  such  stochastic  problems. 

To  improve  the  slow  convergence,  optimal  gPC  bases  are  obtained  by  the  K-L 
expansion  of  the  covariance  matrix  of  the  stochastic  output  solutions  computed  by 
the  adaptive  ANOVA-based  sparse-grid  PCM.  In  the  online  computation,  a  Galerkin- 
projection  based  gPC  method  with  the  optimal  bases  developed  in  the  offline  compu¬ 
tational  part  is  employed,  which  greatly  improves  the  convergence  rate.  The  obtained 
results  in  numerical  examples  considered  indicate  following  three  advantages  of  the 
proposed  adaptive  ANOVA-based  DDSM  method;  (1)  by  integrating  with  adaptive 
ANOVA  decomposition,  it  can  effectively  solve  stochastic  problems  within  desire  ac¬ 
curacy  even  for  problems  with  high-dimensional  and  large  variance  stochastic  inputs; 
(2)  the  same  optimal  bases  can  be  used  for  various  deterministic  forcing  terms  on  the 
right-hand-side  function  of  the  elliptic  PDE  with  random  coefficients;  (3)  comparing  to 
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the  classic  gPC  method,  the  number  of  gPC  modes  to  achieve  specified  accuracy  used 
in  the  adaptive  ANOVA-based  DDSM  method  is  much  smaller,  which  greatly  reduces 
the  computational  cost. 
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